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Summary. We compare different gravitational-wave extraction methods used in three- 
dimensional nonlinear simulations against linear simulations of perturbations of spherical 
spacetimes with matter. We present results from fully general-relativistic simulations of a 
system composed by an oscillating and non-rotating star emitting gravitational radiation. 
Results about the onset of non- linear effects are also shown. 



1 Simulation setup and wave extraction 

We considered a stable equilibrium star model perturbed with an axisymmetric quadrupolar 
fluid perturbation. Assuming the perfect fluid approximation for the matter, the star is 
described by a polytropic equation of state of the typ43 P = Kp'^ with K = 100, F = 2 
and central rest-mass density pc = 1.28 x 10~^. The gravitational mass and the radius are 
M = IAMq and R = 9.57. This model has been widely used in the literature (see e.g. |10)). 
The perturbation consists in an axisymmetric pressure perturbation: the angular pattern is 



We use standard dimensionless units c — G — Mq = 1 unless otherwise specified, for 
clarity. 
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proportional to the (2, 0) spherical harmonic, while the radial profile has no nodes in the star 
interior. Different perturbation amplitudes are considered: A — {0.001, 0.01, 0.05, 0.1}, named 
Ao,. . ., A3. After imposing the fluid perturbation, the linearized (Hamiltonian) constraint is 
solved under the hypothesis of conformal flatness to obtain consistently the perturbed metric. 
See [1] |3] for a detailed description and discussion. 

The same initial data were then evolved with two different codes. The PerBaCCo code, 
described in [3], was used to solve in the time domain the (1+1) equations for the linear 
(even parity) perturbations of a spherical spacetime with matter. The resolution of the 
(radial) ID grid is Ar = 0.032, which corresponds to 300 points inside the star. Gravitational 
waves (GWs) are directly computed via the Zerilli-Moncrief function !f'^''-'(t) (see below). 
The Cactus-Carpet-CCATIE-Whisky code was instead used to solve the complete three- 
dimensional (3D) nonlinear hydrodynamics equations coupled with Einstein's field equations 
within in the "3+1" formalism. See [5] for a general description of the code and the recent 
results obtained. The Cartesian 3D grid consists in 3 fixed refinement levels (boxes) with 
spacing (in each direction) Ax = 0.5 or Ax — 0.25 or Ax = 0.125 for three resolutions used. 

The extraction of gravitational radiation from a numerically computed spacetime is a 
non-trivial task . The current state-of-the-art methods employed by the numerical-relativity 
community are the metric extraction d la Abrahams-Price |4] and the curvature Newman- 
Penrose extraction via the Weyl scalar -04 [S]. If /i^'" are the multipoles of the GW degrees of 
freedom, the metric extraction provides the Zerilli-Moncrief ^'^'^(t) and the Regge- Wheeler 
'^imi^) functions which are directly related to /i^'": 

where A'^ — {£ + 2){£ + !)£{£ — 1). On the other hand if the curvature extraction is used, 
the GW multipoles must be recovered from the Weyl scalar ipi"^(t) using the relation 

/i''"(f) = Vi'"(i) . (2) 

Because of the finite size of the numerical grid, GWs wave are extracted at finite radius. 

In the following we analyze the performances of these methods in the "simple" system 
composed by a single oscillating star. All the results are compared with the ID code which 
gives more accurate data for the waves. Our work, while differing from both, is motivated 
by the same ideas of [71[^, namely testing the GW extraction methods with perturbative 
results and using an oscillating neutron star as a test-bed. It represents the first complete 
analysis of such a kind employing simulations with the Whisky code. 



2 ID waveforms: wave identikit 

We give here a detailed "identikit" of the waves emitted by the system considered and 
computed using the ID perturbative simulations. This waves will serve as a starting point 
to compare with nonlinear simulations described in the following sections. 

The data we are presenting came from a very long simulation with final time of about 1 s. 
The full waveform is plotted in the top panel of Fig. [1] It can be described by two different 
phases: (i) an initial transient, of about half a GW cycle, say up to t ~ 3 ms, related to the 
setup of the initial data (see the black solid line in the left-middle panel), followed by (ii) 
a quasi-harmonic oscillatory phase, where the matter dynamics are reflected into the stellar 
quasi-normal modes (QNMs): 
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Fig. 1 Analysis of the waves produced with the ID perturbative code. The top panel 
shows the entire waveforms !f'2o''(i)- In the middle panels a blow-up of the waveform (black 
solid line) and the fit (blue dashed line) are shown at early (left) and late (right) time: a 
single fit with the QNMs template in Eq. [3] and two modes perfectly reproduces the entire 
waveform. The bottom panels show the power-spectrum density taken, respectively, from 
the waveform's time series in the interval At ~ [0,0.9] s (left, in green) and in the interval 
At ~ [0.9, 1] s (right, in red): the frequency of the pi mode is clearly present only at late 
time. 



'I'^o (t) ~ ^ An cos {2nVnt + (t>n) exp ( — Qfnt) 



(3) 



From the Fourier spectrum of ^'20^(*) ^^^"^ ^ time interval from 3 to about 30 — 90 ms 
(see left-bottom panel), we found that the signal is dominated by the fundamental mode of 
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r/M 

Fig. 2 Analysis of finite-radius extraction effects on ID perturbative waveforms. The left 
panel shows the maximum of ^20' i^) extracted at different radii. The right panels show the 
result of a linear fit with the template in Eq. |4l 

oscillation of the star, the /-mode, at frequency vo = 1581 Hz, with a very weak contribution 
of the first pressure or p-mode (at frequency around vi ~ 3724 Hz). The accuracy of these 
frequencies obtained from Fourier analysis on such long time series has been checked in 
Ref. [3] and is better than 1%. The frequency of the /-mode agrees with that of Ref. [10) 
within 2%. On the interval considered the waves can be perfectly represented by a one-mode 
expansion, i.e. Eq. |3]with A'^ = 1. As a result of a standard least-squares nonlinear fit with 
a one-mode template we obtained the frequency vo ~ 1580.79 ± 0.01 Hz, which is perfectly 
consistent with that obtained via Fourier analysis. For the damping time, we estimated 
ao = 3.984 ± 0.066 sec"^ and thus to = a^^ ^ 0.25 s. The global residual norm) of the 
fit is of the order of 10~®. 

If we consider the entire duration (Is) of the signal it is clear that a one-mode expansion 
is not sufficient to reproduce accurately the waveform (see Figure [T]). The Fourier analysis of 
the whole waveform reveals the presence of the pi-mode, which has longer damping time, and 
the entire signal has been fitted with a two-mode QNM expansion. The estimated values for 
the frequency of the /-mode and the pi-mode are vo = 1583. 737± 0.002 and ui = 3706.941 ± 
0.001, respectively, and the corresponding dumping times tq = 0.268 s and ti = 2.37 s, with 
errors respectively of the order of 0.1% and 2%. The frequencies are consistent with those 
computed via Fourier analysis as well as via the dumping times. This result is particularly 
significant for a time-domain code. 

Figure [T] shows the extraction result at (isotropic) radius of r = 80Af, which is the same 
extraction radius used for the comparison with the 3D results. However in the ID simulations 
we can easily reach up to r = 500M. We checked (see Figure[2]) the convergence of the waves 
with the extraction radius using as a reference point the maximum oi^^o'it). This point can 
be accurately fitted, as a function of the extraction radius, with 

max(^<°)(t,r)) ~a°° + ^ . (4) 

The extrapolated quantity a°° allows an estimate of the error related to the extraction at 
finite distance. The difference in amplitude of a wave extracted at r > 200M with respect to 
a wave extracted at infinity is less than 2% and frequency seems not to be affected (within 
our accuracy) by finite extraction effects. 
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Fig. 3 Comparison between perturbative and 3D waveforms extracted both with the 
Abrahams-Price metric method and Newman-Penrose curvature method. The results refer 
to the evolutions of the initial data with the smallest perturbation amplitude Aq — 0.001. 
The left panel shows ^20 i^) early times, while the right panel shows the entire time series. 
See text for explanations. 

3 3D waveforms: comparison in the linear regime 

We focus first on the data relative to the smallest initial perturbation amplitude A = Aq = 
0.001. The 3D simulations were performed for a simulation time of about 20 ms. Figure |3] 
summaries the results of our findings. The black solid line is the waveform computed with 
the perturbative code and presented in the previous section. We assume it to be the exact 
one. The waveform extracted from the 3D simulation with the Abrahams-Price method is 
the blue dashed line in the figure. The only difference with the perturbative result is the 
presence of a "burst" of radiation at the early times (see left panel) 0.5 < t < 1 ms. In pP 
we discussed in detail this feature of the metric wave extraction. The main conclusions are 
the following, ('ij It is unphysical: it linearly grows with the extraction radius and produces 
a violation of the linearized Hamiltonian constraint at the wave-extraction level, (ii) It is 
related to an inaccurate implementation of the Zerilli-Moncrief function in the 3D code and 
it can be reproduced in the perturbative code introducing "by hand" inaccuracies (very 
low resolution, low-order approximations for finite differences), (in) It is absent in vacuum 
simulations, e.g. binary black holes, where higher grid resolutions can be used. Apart from 
the junk the metric waveform is perfectly superposed (until f ~ 7 ms) to the exact one; 
frequency and amplitude are consistent within the error bars. 

The GWs extracted via the ^4 scalar require the solution of Eq. [2]to obtain the GW mul- 
tipole (or equivalently the Zerilli-Moncrief function) . The procedure we followed is described 
in Ref. [TT] and consists in computing a double integral of ipi"^ (the "bare" metric waveform) 
and then correcting it by subtracting a polynomial floor. The red dashed-dotted line in Fig. |3] 
is the "bare" waveform; as suggested by the figure, and differently from Ref. [TT], we found 
that a quadratic polynomial is here necessary to align the waveform (see the green dashed- 
dotted line). Once this tuning has been done the -04 waveform is exactly superposed to the 
perturbative one and, differently from the metric waveform, it is free from the unphysical 
burst of junk radiation. 

In summary we can conclude that, in the linear regime, both extraction methods give 
a very good agreement with the "exact" waveform but neither is free from drawbacks. On 
one hand the metric extraction contains a burst of junk radiation at early times that can 
compromise the study of phenomena in the first 5 ms of simulation (e.g. the neutron-star 
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Fig. 4 Nonhnear effects in the waves computed with the 3D code. The left panel shows 
the amplitude of the (2,0) mode, determined by fitting Eq. [3]with N — 1, as a function of 
the initial perturbation amplitude A; while the linear scaling of the linear code is perfect, in 
the nonlinear case the deviations, that can be seen already for small values of A = Ai = 0.01, 
become conspicuous for A > Ai. The right panel shows the power-spectrum density of the 
integrals of the rest-mass density projected trough the (2, 0) spherical harmonic. For larger 
values of A more frequencies are present because of nonlinear mode couplings. 

collapse to black hole or neutron-star ui-modes) . On the other hand the curvature extraction 
requires tuning of the integration constants. 



4 3D waveforms: the nonlinear regime 

Finally we discuss nonlinear effects in the 3D waveforms. Previous studies can be found for 
example in [S] (second-order perturbation theory) and in [S] (CFG approximation of Einstein 
equations). We confirm those results without using approximations. Referring to the left 
panel of Fig. [J] we found that the amplitude of the (2, 0) mode (computed fitting Eq. ^ 
scales linearly with the initial perturbation amplitude A - as expected -, but it progressively 
decreases in the nonlinear case due to excitation of other channels. The nonlinear effects are 
already visible for A = Ai (10% difference in A20) and quickly increase for A2 (38%) and 
A3 (66%). We note that the frequency of the wave, instead, is essentially not affected by 
nonlinearities with our setting. 

Furthermore we considered the Fourier transform of the rest-mass density projected 
(multiplied and integrated on the 3D-space) on the spherical harmonic (2,0). As the right 
panel of Fig. |4] shows, for A = Ao (green line) only the £ — 2 f and pi (and P2) mode 
frequencies are present, while for higher values of A the spectra are rich of other frequencies. 
Most of the frequencies can be associated to the radial {£ = 0) fundamental frequency F 
(1466 Hz), its overtone Hi (5955 Hz), and to combination tones between the radial and the 
£ — 2 frequencies: ^'nniin = ^0 ^^2- This analysis prove that mode couplings are present in 
the dynamics of the system already for small initial perturbation amplitudes, but they affect 
strongly only the amplitude of the GWs. 
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